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ABSTRACT 

If the accelerated expansion of the Universe at the present epoch is driven by a dark en- 
ergy scalar field, there may well be a non-trivial coupling between the dark energy and the cold 
dark matter (CDM) fluid. Such interactions give rise to new features in cosmological structure 
growth, like an additional long-range attractive force between CDM particles, or variations of 
the dark matter particle mass with time. We have implemented these effects in the N-body 
code GADGET-2 and present results of a series of high-resolution N-body simulations where 
the dark energy component is directly interacting with the cold dark matter. As a consequence 
of the new physics, CDM and baryon distributions evolve differently both in the linear and 
in the nonlinear regime of structure formation. Already on large scales a linear bias develops 
between these two components, which is further enhanced by the nonlinear evolution. We also 
find, in contrast with previous work, that the density profiles of CDM halos are less concen- 
trated in coupled dark energy cosmologies compared with ACDM, and that this feature does 
not depend on the initial conditions setup, but is a specific consequence of the extra physics 
induced by the coupling. Also, the baryon fraction in halos in the coupled models is signifi- 
cantly reduced below the universal baryon fraction. These features alleviate tensions between 
observations and the ACDM model on small scales. Our methodology is ideally suited to ex- 
plore the predictions of coupled dark energy models in the fully non-linear regime, which can 
provide powerful constraints for the viable parameter space of such scenarios. 
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1 INTRODUCTION 

The last decade has seen an astonishing amount of new cosmo- 
logical data from many different experimen ts, ranging from large- 
scale structure surveys jPercival et al J 1 2001 , e.g.) to Cosmic Mi- 
crowave Ba ckground (CMB) jKomatsu et al.l 20081) and Type la 
Super novae teiess et al]| 19981 : IPerlmutter et afll 19991 : IXstier et all 
2006) observations. These experiments all consistently show that 
the Universe is almost spatially flat with a current expansion rate of 
about 70 km s~ 1 Mpc~ 1 , and contains a total amount of matter that 
accounts for only ~ 24% of the total energy density. From a theore- 
tical point of view, understanding the remaining 76% - which must 
be in the form of some dark energy (DE) component able to drive 
an accelerated expansion - is a serious challenge. 

A simple cosmological constant would be in agreement with 
a large number of observational datasets, but would also raise two 
fundamental questions concerning its fine-tuned value and the co- 
incidence of its domination over cold dark matter (CDM) only at 
a relatively recent cosmological epoch. An alternative consists of 
identifying the dark energy component with a dynamic scalar field 
dWetterichll 988: Ratr a & Peebles! 19881) . thereby seeking an expla- 



nation of the fundamental problems challenging the cosmological 
constant in the properties of the dynamic evolution of such scalar 
field. 

An interesting suggestion recently developed concerns the in- 
vestigation of possible interactions between the d ark energy scalar 
field and other matter species in the Universe dWettericbl 1 19951 : 
Amendol j|200d: iFarrar & Peeblej|2004l : iGubser & Peeblesil2004l : 
Farrar & RoserJl2007l) . The existence of such a coupling could pro- 
vide a unique handle for a deeper understand i ng of the DE problem 



dOuartin et al.l 



Manga no et al 



20081: [Brookfield et al. 2008]; iGromov et al.1 [20041: 



20031 : lAnderson & Carroll 1997b- It is then cru- 



cial to understand in detail what effects such a coupling imprints 
on observa ble features like, for example, the CMB and structure 
formation ( La Vacca et "aHl20q9l;lBean et al.ll2008 | :lBertolami et all 
| 2007l: iMatarrese et alJ 120031: IWang et alj 120071 : Iguo et alj I2007F 
Main ini & Bonomettol2o'o73 : iLee et al.ll2006h . 

In this paper, we perform the first fully self-consistent high- 
resolution hydrodynamic N-body simulations of cosmic structure 
formation for a selected family of coupled DE cosmologies. The 
interaction between DE and CDM is expected to imprint characte- 
ristic features in linear and nonlinear structures, and could possi- 
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bly open up new ways to overcome a series of observational chal- 
lenges for the ACDM concordanc e cosmology, rangin g from the 
satellite abundance in CDM halos jNavarro et al.ll 1996th to the ob- 
served low baryon fraction in large galaxy clus ters ( I Allen et all 
|2004 IVikhlinin etai]|2006l ; iLaRoque et aitl2006h . and to the so- 
called "cusp-core" problem for the density profiles of CDM halos, 
that was first raised in relat ion to observations of dwarf and low- 
surface-brightness galaxies dMoore|[l994l ; lFlores & Primacldl 19941 ; 
ISimon et al.ll2003l) . but that more recent studies seem to extend to 
a wider range of objects including Milky- Way type spiral galaxies 
(N avarro & Steinmetz 2000; Binney & Evans 2001) up to galaxy 
groups and clusters dSand et alj2002ll2004l ; lNewman et alj[2 009). 

In the following we will consider DE as the energy ass ociated 
with a scalar field (i> Iwetterichll988l : rRatra & PeeblesTl988h whose 
energy density and pressure are defined respectively as the (0, 0) 
and (0, i) components of its stress energy tensor, so fha{]J 



P<t> 



2 a 2 



U{<t>) , 



(1) 



(2) 



with U((j>) being the self interaction potential. The models we in- 
vestigate in this work by means of detailed high-resolution N-body 
simulations are ph ysically identical to the ones already studied by 
Maccioet "al]( l2004h . With respect to this previous work, however, 
our analysis significantly improves on the statistics of the number 
of cosmic structures analyzed in each of our different cosmologi- 
cal models and on the dynamic range of the simulations. We also 
extend the analysis to a wider range of observable effects arising 
from the new physics of the DE-CDM interaction, and we include 
hydrodynamics. Despite the physical models investigated be iden- 
tical, our outcomes are starkly different from the ones found in 
lMacci6etai] ( l2004l) . 

This work is organized as follows. In Section|2]we detail the 
coupled DE cosmologies under investigation, both regarding the 
background evolution and the perturbations, and we introduce a 
numerical package to integrate the full set of background and per- 
turbation equations up to linear order. In Section|3] we describe the 
numerical methods we use, and the particular set of simulations we 
have run. We then present and discuss our results for N-body sim- 
ulations for a series of coupled DE models in Section|4] Finally, in 
Section|5]we draw our conclusions. 



2 COUPLED DARK ENERGY COSMOLOGIES 

Cou pled cosmologies can b e described following the considera- 
tion dKodama & Sasaki 1984) that in any multicomponent system, 
though the total stress energy tensor T M „ is conserved 



(3) 



the T M „( Q ) for each species a is, in general, not conserved and 
its divergence has a source term Q( Q ) M representing the possibility 
that species are coupled: 



;(q)m i 



(4) 



1 In this work we always denote with a prime the derivative with respect 
to the conformal time t, and with an overdot the derivative with respect 
to the cosmic time t. The two time variables are related via the equation: 
dr = dt/a. 



with the constraint 



(5) 



Furthermore, we will assume a flat Friedmann-Robertson- Walker 
(FRW) cosmology, in which the line element can be written as 

ds 2 = a 2 (r)(— dr 2 + 8ijdx z Ax i ) where o(r) is the scale factor. 

2.1 Background 

The Lagrangian of the system is of the form: 
1. 



C = --d^d^ - U(<f>) - m(0)VV + Ad* 



(6) 



in which the mass of matter fields ip coupled to the DE is a function 
of the scalar field (f>. In the following we will consider the case in 
which the DE is only coupled to cold dark matter (CDM, hereafter 
denoted with a subscript c). The choice m((f>) specifies the coupling 
and as a consequence the source term Q(^,) M via the expression: 

d\xim(4>) 



Q 



(0V 



-pa <9 M </>. 



(7) 



Due to the constraint (O, if no other species is involved in the cou- 
pling, Q {chl = -Q W)M . 

The zero-component of equation © provides the conservation 
equations for the energy densities of each species: 

p'^ = —3T-Lp^(l + w^) — Q(</,)o , (8) 
Pc — —3Hp c + Q(4>)o ■ 

Here we have treated each component as a fluid with T"^^ = 
(p a + Pct)Uftu" + p a S^, where u M = (—a, 0, 0, 0) is the fluid 4- 
velocity and w a = p a / p a is the equation of state. The class of 
models considered here corresponds to the choice: 



m((f>) = moe"' 3 ^ , 

with the coupling term equal to 



AI 



it 

-Pc<p ■ 



(9) 



(10) 



This set of cosmologies has been widely investigated, for P(<j>) 
given by a co nstant, both in i t s background and linear perturba- 
tion features dWetterichlll99l : lAmendolal l2000l) as well as with 
regard to the effects on stru cture formation (Amendola 120041 ; 
IPettorino & Baccigalupill2008l) , and via a first N-body simulation 
dMaccio et al. 2004). 



2.2 Linear perturbations 

We now perturb the quantities involved in our c osmological frame- 
work up to the first order in th e perturbations (Kod ama & Sasakj 
ll984l ; lMa& Bertschinger 1995). The perturbed metric tensor can 
then be written as 



<7 M „(T,X) = g AW (T)+<5g A .„(T,x) 



(ID 



where Sg^ <C 1 is the linear metric perturbation, whose expres- 
sion in Fourier space is given by: 

Sgoo = -2a 2 AY , 8g 0l = -VSY, , 



5 gij = 2a 2 [H L Y8 ij + H T Y id 



(12) 



where A, B, Hl, Ht are functions of time and of the wave vec- 
tor k, and Yi , Yij are the vector and tensor harmonic functions ob- 
tained by deriving Y, defined as the solution of the Laplace equation 
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S'i'Vi'VjY = — |fc| 2 y. Analogously, the perturbed stress energy 
tensor for each fluid (a) can be written as T, M , = T? , + <5T, M > 
where the perturbations read as: 

ST ( a )0 =-P(a)S( a )Y, (13) 

= V) («(..) --B)^, 

5T (a)0 = ~h('*) V (°') Y ' > 

The perturbed conservation equations then become: 

(p a S a )' ' +3Hp a S a +h a (kv a +3H' L )+3'Hp a TTL a = -8Q( a ) (14) 
for the energy density perturbation 8 a — 8p a j p a , and: 



[h a {v a - B)\' + mh a (v a - B) 
2 

-kp a TlLa — fl a kA + -kpiVTa = <5<3(a)i 



(15) 



for the velocity perturbation u a . 

The scalar field <f> can also be perturbed, yielding in Fourier 
space 

4> = <f> + 5(f> = <(> + X (t)Y . (16) 
Furthermore, we can express the perturbations of the source as: 

/9(0 -,, P(4>),,, fo., 



8Q(<j,)a — - 
5Q{4>)i = k 



M 
M 



■p c S(f> 



(17) 



(18) 



In the Newtonian gauge (B = 0, H T — 0, H L — *, A = <J>) the 
set of equations for the density and velocity perturbations for DE 
and CDM read: 



8p'^ + 3 / H{8p 4> + 8pt) + kh^vj, + 3/ty*' = (19) 

~JF P ^ ~~M~^ Pc ~M~ ^ Pc ' 
5p' c + 3H5p c + kp c v c + 3p c & = 

~^r pM ~ ~m~^ Spc ~ ~~m Hpc > 

h^v'^ + (h'^ + 4Uh<t,) V4, — k8p<f, — kh^ = k ^^- pcScf) , 



H ^Ls' ]v c -k* = -k 



M 



M r J M 
The perturbed Klein Gordon equation in Newtonian gauge reads: 

5<t>" + + (k 2 + aU^) 5<t> - </>' (V - 3*') + 



2a [7 * = 3H 2 fl c [J3(<j>)5c + 2/3(0)* + (3^ 



(20) 



For the N-body implementation we are interested in, the New- 
tonian limit holds, for which A = H/k <C 1. In this case we have 

8(f> ~ 3\ 2 fl c f3{(t))S c . (21) 

In this limit, the gravitational potential is approximately given by 

, 3 A 2 



2M 2 



We can then define an effective gravitational potential 



<£ c = * + 



M 



(22) 



(23) 



which also reads, in real space and after substituting the expressions 
for * (Eqn.HUl and for 8(f) (Eqn.|2T): 



V 2 *c = - \po8 c (1 + 2/3 2 (0)) - y ^ p a 8 a , 



(24) 



where the last term takes into account the case in which other com- 
ponents not coupled to the DE are present in the total energy budget 
of the Universe. Cold dark matter then feels an effective gravita- 
tional constant 



G c = Giv[l + 2/3 2 (0)], 



(25) 



where Gjv is the usual Newtonian value. Therefore, the strength 
of the gravitational interaction is not a constant anymore if /3 is a 
function of the scalar field (f>. The last equation in l !19t , written in 
real space and in terms of the effective gravitational potential, gives 
a modified Euler equation of the form: 

Vt) c + rt —<P I V1)c 



M 



n c 8 c + 2n c 8 c /3 2 (<t>) + 



= 0. 



(26) 



As in lAmendolj d2004h . if we assume that the CDM is con- 
centrated in one particle of mass m c at a distance r from a particle 
of mass M c at the origin, we can rewrite the CDM density contri- 
bution as 



£l C 8r = 



8nGM c e 



(27) 



3ft 2 a 

where we have used the fact that a non-relativistic particle at po- 
sition r has a density given by m c nS(r) (where S(r) stands for 

the Dirac distribution) with mass given by m c oc e J 8 ^' d '* j for- 
mally obtained from equation l(9](. We have further assumed that the 
density of the M c mass particle is much larger than p c . The Euler 
equation in cosmic time (dt — a dr) can then be rewritten in the 
form of an acceleration equation for the particle at position r: 



= -Hv c - V 



G C M C 



(28) 



where we explicitely see that the usual equation is modified in three 
ways. 

First, the velocity-dependent term now contains an additional 
contribution given by the second term of the expression defining 
H: 



H = H 1 



M H 



(29) 



Second, the CDM particles feel an effective gravitational constant 
G given by J25t . Third, the CDM particles have an effective mass, 
varying with time, given by: 



AL 



Mc e-I^ da 



(30) 



Eqn.|28]is very important for our discussion since it represents the 
starting point for the implementation of coupled DE models in an 
N-body code. We will discuss in detail how this implementation is 
realized in Sec. 13. II but it is important to stress here that Eqn.1281 
is written in a form that explicitely highlights its vectorial nature, 
which has not been presented in previous literature. The vectorial 
nature of Eqn. [28] is a key point in its numerical implementation 
and therefore needs to be properly taken into account. 
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In the N-body analysis carried out in the present paper, we 
consider /3 to be a constant, so that the effective mass formally reads 
M c = M c e _/3 (* _ *o) w e defer an investigation of variable /3(</>) to 
future work. We have numerically solved the full background and 
linear pertu rbation equ ations with a suitably modified version of 
CMBEASY ( Doran 2005), briefly described in the following section. 



2.3 Background and linear perturbation integration of 
coupled DE models: a modified CMBEASY package 

We have implemented the full background and linear perturbation 
equations d erived in th e previous section in the Boltzmann code 
CMBEASY (Doran 2005) for the general case of a DE component 
coupled to dark matter via a coupling term given by Eqn. {4}. The 
form of this coupling, as well as the evolution of the DE (either 
modelled as a scalar field or purely as a DE fluid), can be freely 
specified in our implementation. 

Compared to the standard case of uncoupled DE, the modifi- 
cations include a modified behaviour of the background evolution 
of CDM and DE given by Eqs. ([8j and ((9}, as well as the implemen- 
tation of the linear perturbations described by Eqn. H9I >, and their 
corresponding adiabatic initial conditions. The presence of a CDM- 
DE coupling furthermore complicates the choice of suitable initial 
conditions even for the background quantities, since dark matter 
no longer scales as a~ 3 , and so cannot simply be rescaled from its 
desired value today. For each of the models considered here, we 
choose to set the initial value of the scalar field close to its tracker 
value in the uncoupled case, and then adjust the value of the poten- 
tial constant A (see Eqn. \3ll below) and the initial CDM energy 
density such that we obtain the desired present-day values. 



Parameter 


Value 


^CDM 


0.213 


H 


71.9 km s" 1 Mpc -1 




0.743 


C8 


0.769 


fib 


0.044 


n 


0.963 



Table 1. Cosmological parameters for our set of N-body simulations, 
consistent with the WMAP 5 year results for a ACDM cosmology 
iKomatsu et alj2008t) . 

for a ACDM cosmological model. In the coupled models we con- 
sider, the role of DE is played by a qu intessence scalar field with a 
Ratra-Peebles dRatra & Peeb les 1988) self-interaction potential of 
the form: 

A 4+Q 



where A and a fix the DE scale in agreement with observations, 
and with a constant coupling to CDM particles only, a s described in 
Sec. | 2 . 1 1 we label them as RP1-RP6 in analogy with lMaccio et al.l 
d2004h . 

For four of these models (ACDM, RP1, RP2, RP5) we then ran 
high-resolution simulations in a smaller cosmological box (Lbox = 
80ft -1 Mpc, iVpart = 2 x 512 3 ), and we investigated the proper- 
ties of collapsed objects for this set of simulations. In addition to 
these four high-resolution simulations we ran other three simula- 
tions with the same resolution (ACDM-NO-SPH, RP5-NO-SPH, 
RP5-NO-GF), whose features will be described below. The cos- 
mological parameters for our models are listed in Table Q] and the 
physical parameters of each model together with the details of the 
corresponding N-body simulations are listed in Table[2] 



3 THE SIMULATf ONS 

The aim of this work is to investigate the effects that a coupling 
between DE and other cosmological components, as introduced in 
Sec. [2] can have on cosmic structure formation, with a particular 
focus on the nonlinear regime which is not readily accessible by 
the linear analytic approach described in Sec. 12.21 To this end we 
study a set of cos mological N-bo dy simulations performed with the 
code GADGET-2 ( Springel 2005), that we suitably modified for this 
purpose. 

The required modifications of the standard algorithms of 
the N-body code for simulating coupled DE cosmologies are 
extensively described in Sec. 13.11 Interestingly, they require to 
reconsider and in part drop several assumptions and approx- 
imations that are usually adopted in N-body simulations. We 
note that previous attempts to use cosmological N-body simula- 
tions for different flavours of modified Newtonian gravity have 
been d i scussed, for example, in [Macci6 et al. (2004); Nusse ret al.l 
j2005l):[stabenau & Jainl j200o1):lKesden & Kamionkowskil j2006h: 
ISpringel & Farraj feOfM); iLaszlo & Beanl J2008I) : ISutter & Ricked 
(2008)jlOyaiztj fe008h: llCeselman et alj d2009r) but to our knowl- 



edge iMaccio et alj d2004h is the only work to date focusing 
on the properties of nonlinear structures in models of coupled 
quintessence. 

Therefore, with our modified version of GADGET-2 we first 
ran a set of low-resolution cosmological simulations (Lbox = 
320/i" 1 Mpc, jV part = 2 x 128 3 ) for the same coupled DE models 
investigated in M accio et a I]( l2004h . but with cosmolo gical parame- 
ters updated to the latest results from WMAP dKomatsu et al.l2 008) 



3.f Methods 

The presence of a direct coupling between the DE scalar field (j> 
and other cosmic fluids - in the fashion described by Eqs. {4] [8] 
[9]l - introduces new features in the cosmic evolution as well as 
additional physical processes that must be taken into account in 
N-body models. In the following, we describe these features and 
their implementation in GADGET-2 one by one, recalling an d fur- 
th er emphasising the results descr ibed in lMaccio" et al.l J2004h and 
m lPettorino & Baccigalupl d2008h . 



3.1.1 Modified expansion rate 



As pointed out in Sec. 12.11 the coupling modifies the background 
evolution through the e xistence of a phas e - corresponding to the 
so called <^>MDE era in Amendola (2000) - in which the coupled 
matter fluid (CDM in our case) and the DE scalar field evolve with 
a constant energy density ratio (here we always assume the Uni- 
verse to be flat such that fi t ot = 1). This lea ds to the presence 
of a non negligible Ear ly DE component (EDE. lDoran et al]|200ll : 
iDoran & "Ro bbers 2006) during the entire epoch of structure forma- 
tion. The effect of such an EDE is to change the expansion history 
of the Universe, which has to be properly taken into account for the 
N-body time integration. In order to do so, we replaced the com- 
putation of the Hubble function H(a) in GADGET-2 by a linear 
interpolation from a table of values of H(a) precomputed for each 
model under investigation with the modified version of CMBEASY 
described above. The effect of the coupling on the expansion is 
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Model 


a 


fa 


fa 


Box Size (ft -1 Mpc) 


Number of particles 


M b (h- 1 M Q ) 


Modm (h~ l M Q ) 


e s (ft, -1 kpc) 


ACDM (low) 











320 


2 X 128 a 


1.9 x 10 11 


9.2 x 10 11 


50.0 


ACDM (high) 











80 


2 X 512 3 


4.7 X 10 7 


2.3 x 10 s 


3.5 


ACDM (high - no SPH) 











80 


2 X 512 3 


4.7 X 10 7 


2.3 x 10 s 


3.5 


RP1 (low) 


0.143 





0.04 


320 


2 X 128 3 


1.9 X 10 11 


9.2 x 10 11 


50.0 


RP1 (high) 


0.143 





0.04 


80 


2 X 512 3 


4.7 X 10 7 


2.3 x 10 s 


3.5 


RP2 (low) 


0.143 





0.08 


320 


2 X 128 3 


1.9 x 10 11 


9.2 x 10 11 


50.0 


RP2 (high) 


0.143 





0.08 


80 


2 X 512 3 


4.7 x 10 7 


2.3 x 10 s 


3.5 


RP3 (low) 


0.143 





0.12 


320 


2 X 128 3 


1.9 x 10 11 


9.2 x 10 11 


50.0 


RP4 (low) 


0.143 





0.16 


320 


2 X 128 3 


1.9 x 10 11 


9.2 x 10 11 


50.0 


RP5 (low) 


0.143 





0.2 


320 


2 X 128 3 


1.9 x 10 11 


9.2 x 10 11 


50.0 


RP5 (high) 


0.143 





0.2 


80 


2 X 512 3 


4.7 x 10 7 


2.3 x 10 s 


3.5 


RP5 (high - no SPH) 


0.143 





0.2 


80 


2 X 512 3 


4.7 X 10 7 


2.3 x 10 8 


3.5 


RP5 (high - no GF) 


0.143 





0.2 


80 


2 X 512 3 


4.7 X 10 7 


2.3 x 10 s 


3.5 


RP6 (low) 


2.0 





0.12 


320 


2 x 128 3 


1.9 x 10 11 


9.2 x 10 11 


50.0 



Table 2. List of the different simulations performed with our modified version of GADGET-2. The simulations have different force and mass r esolution, and 
are ac cordingly labelled as low or high resolution. Notice that the cited values of the coupling listed here are different from the ones adopted in Maccio et al. 
12004 ) due to the d ifferent definition of the coupling function (9)- However, the models in effect have identical coupling strength to those investigated in 
iMaccio etai]fe004l) . 



Function 


Meaning 


H(a) 


Hubble function 


AG(a) 


Possible global variation of the gravitational constant 


fa{4>) 


Coupling function for the baryons 


fa{4>) 


Coupling function for CDM 


Am b 


Correction term for baryon particle masses 


Am c 


Correction term for CDM particle masses 




Dimensionless kinetic energy density of the scalar field 



Table 3. List of input functions for the coupled DE implementation in 
GADGET-2. 



shown in Fig. Q] We note that the same approach has also been 
adopted for the other relevant quantities described in Table[3] which 
were first computed numerically using CMBEASY, and then used as 
an input for our modified version of GADGET-2. 



3.1.2 Mass variation 

As described in Sec. 12.11 the coupled species feature an effective 
particle mass that changes in time. Consequently, the correspond- 
ing cosmological densities p c or p b will not be scaling anymore 
as a -3 , but will have a correction term arising from the variation 
of particle masses on top of the pure volume dilution. This correc- 
tion depends on the scalar field dynamics, and takes the form of 
Eqn. d30t . Then, if the particle masses in the simulation box are 
computed according to the cosmic densities fl Ci o and Q.b,o at the 
present time, they will need to be corrected at each timestep by a 
factor 

Am c (a) = e Ja do , (32) 

Am b (a) = e -J> w ^ d \ (33) 

In Fig. [2] we show the evolution with redshift of the correction 
term d32b for our set of models. 

Although the coupled DE implementation for GADGET-2 that 
we present here in principle allows for a coupling also to baryons, 
such a coupling is tightly cons trained by experime ntal tests of 
gravity on solar system scales toamoureta"uT l990). Therefore, 
in the coupled DE models we consider in this work, no cou- 
pling to baryons is considered, and baryon masses are always 
constant. However, even a very tiny coupling to baryons, in the 



Mass correction for different coupled dark energy models 



1.4 



1.0 





ACDM 






RP1 






RP2 






RP3 






RP4 






RP5 






RP6 
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Figure 2. Mass correction as a function of redshift for the coupled DE mod- 
els with constant coupling. 

range allowed by present bounds, could possibly play a significant 
role in cosmological m odels with multiple dark matter families 
fauev&Wandel3l2006l) . like for example the Growing Neutrino 
scenario introduced by lAmendola et alTd2008l) . We plan to study 
detailed N-body simulations of this kind of models in forthcoming 
work. 



3.1.3 Cosmological extra velocity-dependent acceleration 

A further modification to the GADGET-2 code concerns the extra 
cosmological velocity-dependent term induced by the coupling in 
the Euler equation, as shown in Equations l |28l \29\ . In standard 
cosmological simulations with GADGET-2 the usual cosmological 
velocity-dependent term (first term of Eqn.|29l> is not directly com- 
puted because the choice of variables removes it from the accelera- 
tion equation. If we denote with x comoving coordinates and with 
r — a(t)x physical coordinates, we have that: 

r — Hr + v p , v p = a(t)x . (34) 

Instead of using the peculiar velocity v v , GA DGET-2 makes use of 
the variable p = a 2 (t)x (see lSpringell J2005I) for more details), for 
which the following relation holds: 
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Figure 1. Left panel: Hubble functions as a function of redshift for the different coupled DE models with constant coupling investigated in this work and 
described in Table[2]as compared to ACDM (black curve). Right panel: Ratio of the Hubble functions of each coupled DE model to the Hubble function for 
the reference ACDM cosmology as a function of redshift. 



1 . H 
v p = -p p . 

a a 



(35) 



It is then straightforward, by using Eqn. d35l l. to find a generaliza- 
tion of Eqn. d28t to a system of N particles, in terms of the new 
velocity variable p: 



1 



\ ~< GijlTbjXi 



(36) 



where i and j are indices that span over all the particles of the 
simulation, 7 = c, 6 for CDM or baryons respectively, and Gij 
is the effective gravitational constant between the i-th and the j- 
th particles, as determined in Eqn. ( 1251 ) and whose implementation 
will be discussed below. 

It is evident from Eqn. \36\ that for zero coupling no cosmo- 
logical velocity-dependent term is present in the acceleration equa- 
tion, which is then just Newton's law in comoving coordinates. In 
general, however, whenever a coupling is present, the additional 
term 



a(p + Ap) = ap + P 1 {cj))—ap i 



(37) 



has to be explicitely added to the Newtonian acceleration of every 
particle. This term does not depend on the matter distribution in the 
rest of the Universe. It is therefore a purely cosmological drag that 
would be present also in absence of any gravitational attraction. It 
is interesting to notice that in the case of constant positive coupling 
and a monotonic scalar field potential as investigated here, the ex- 
tra cosmological velocity-dependent term induces an acceleration 
in the same direction as the velocity of the particle. It therefore is 
effectively a "dragging" term, speeding up the motion of the parti- 
cles. 

Scalar field models where the dynamics of the field or the 
evolution of the coupling induce a change in the sign of /3i(<p)<fi, 
thereby changing the direction of th is extra velocity-depe ndent 
force, will be studied in future work tealdi & Mac cio in prep.). 



3.1.4 Fifth force implementation 

One of the most important modifications introduced by the cou- 
pling between DE and CDM is the presence of a modified gravi- 



tational constant, formally written as in Eqn. ( 1251 ), for the gravita- 
tional i nteraction of CD M particles. In fact, if in general the substi- 
tution ( Amendola 200jl) 



Gn — > Gim = Gn ■ (1 + 2/3;/3 m ) . 



(38) 



holds for each pair (I, m) of particles, with I and m denoting the 
species of the particle, in our case only CDM-CDM interaction is 
affected, while baryon-CDM or baryon-baryon interactions remain 
unchanged since fib = 0. The dependence of this modified gravi- 
tational interaction on the particle type requires an N-body code to 
distinguish among different particle types in the gravitational force 
calculation. In GADGET-2, the gravitational interaction is co mputed 
by means of a TreePM hybrid method (see Springel 2005, for de- 
tails about the TreePM algorithm), so that both the tree and the 
particle-mesh algorithms have to be modified in order to account 
for this new phenomenology. 

Tree algorithm modifications - In a standard tree algorithm, 
each node of the tree carries informations about the position of its 
centre of mass, its velocity, and its total mass. The decision whether 
to compute the force exerted on a target particle by the whole node 
or to further divide it into eight smaller nodes is made based on a 
specific opening criterion, which sets the accuracy threshold for ap- 
proximating the gravitational potential of a distribution of particles 
with its low-order multipole expansion. Since in uncoupled cos- 
mological models all particles interact with the same gravitational 
strength, as soon as the opening criterion is fulfilled the force is 
computed assigning all the mass contained in the node to its centre 
of mass. For coupled quintessence cosmologies, this is no longer 
accurate enough given that the different particle species will con- 
tribute differently to the gravitational force acting on a target par- 
ticle. This means that besides the total mass and the total centre 
of mass position and velocity, each node has to carry information 
about the mass and centre-of-mass position and velocity of each 
particle species with different coupling. 

Particle-Mesh algorithm modifications - In the Tree-PM al- 
gorithm, the long-range part of the gravitational force is computed 
by means of Fourier techniques. For coupled DE models, where 
different particle species interact with an effectively different grav- 
itational force, the PM procedure has to be repeated as many times 
as there are differently interacting particle types, each time assign- 
ing to the cartesian grid only the mass in particles of a single type, 
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Figure 3. Matter power spectra at z = 60 for interacting DE models with 
constant coupling as computed by CMBEASY. 



and then computing the gravitational potential and the acceleration 
deriving from the spatial distribution of that particle species alone. 
In this way, the total force is built up as a sum of several partial 
force fields from each particle type. 



3.1.5 Initial conditions 

The initial conditions of a cosmological N-body simulation need to 
specify the positions and velocities of all the particles in the cos- 
mological box at the starting redshift z, of the simulation. These 
quantities are usually computed by setting up a random-phase re- 
alization of the power spectrum of the studi ed cosmological m odel 
according to the Zel'dovich approximation t Zel'dovichl 1970 ). The 
normalization amplitude of the power spectrum is adjusted such 
that the linearly evolved rms-fiuctuations erg on a top-hat scale of 
8 ft -1 Mpc at a given redshift z„ rm (usually chosen to be z norm = 
0) have a prescribed amplitude. 

The coupling between DE and CDM can have a strong im- 
pact on the transf er function of matter density fluctuations, as 
first pointed out bv lMainini & Bonomettol J2007bh . For this reason 
we compute the required initial power spectrum directly with the 
modified Boltzmann code CMBEASY, because the phenomenolog- 
ical parameterizations of th e matter power spectrum available for 
the A CDM cosmology (e.g. iBardeen et~ai]|l986l : lEisenstein & Hvj 
1998) would not be accurate enough. The resulting effect on the 
power spectrum is shown in Fig.[3]for the different models consid- 
ered in our set of simulations. 

Once the desired density field has been realized with this pro- 
cedure, the displacements of the particles from the grid points need 
to be rescaled with the linear growth factor D+ for the cosmolog- 
ical model under investigation between the redshifts z norm and Zi 
in order to set the correct amplitude of the power spectrum at the 
starting redshift of the simulation. Also, the velocities of the par- 
ticles are related to the local overdensities according to linear per- 
turbation theory, via the following relation, here written in Fourier 
space: 



v(k, a) = if(a)aH8(k, a] 



where the growth rate /(a) is defined as 
d\nD+ 



/(«) 



din a 



(39) 



(40) 



This requires an accurate calculation of the linear growth function 
D+(z) for the coupled model, which we again compute numeri- 
cally with CMBEASY. 

We note that a phenomenological parameterization of the 
growth function for coupled DE models with constant cou- 
pling to dark matte r has recently been made available by 
|Pi Porto & Amendolal j2008l) . However, it is only valid for mod- 
els with no admixture of uncoupled matter, whereas in our case we 
also have a baryonic component. Also, in the ACDM cosmology, 
the total growth rate is well approximated by a power of the total 
matter density fijp with 7 = 0.55, roughly independently of the 
cosmological constant density jPeebleslI 1 980h . This is however no 
longer true in coupled cosmologies, as we show in Fig. [4] We find 
that, for our set of coupled DE models, a different phenomenologi- 
cal fit given by 



/(a) "02,(1 +7 



' a 2 



(41) 



with 7 = 0.56 (as previously found in I Amendola & Ouercellini 
l2004h and e c = 2.4 works well. The fit gTJi reproduces the growth 
rate with a maximum error of ~ 2% over a range of coupling values 
between and 0.2 and for a cosmic baryon fraction Qb/Q m at 
z — in the interval 0.0 — 0.1 for the case of the potential slope 
a = 0.143 (corresponding to the slope of the RP1-RP5 models). 
For a value of a — 2.0 (corresponding to the slope assumed for 
RP6) the maximum error increases to ~ 4% in the same range of 
coupling and baryon fraction. In Fig. [4] we plot both the fitting 
formulas together with the exact f(a). For our initial conditions 
setup we in any case prefer to use the exact function f(a) directly 
computed for each model with CMBEASY, rather than any of the 
phenomenological approximations. 



3.2 Tests of the numerical implementation: the linear growth 
factor 

As a first test of our implementation we check whether the linear 
growth of density fluctuations in the simulations is in agreement 
with the linear theory prediction for each coupled DE model un- 
der investigation. To do so, we compute the growth factor from 
the simulation outputs of the low-resolution simulations described 
in Table [2] by evaluating the change in the amplitude of the mat- 
ter power spectrum on very large scales, and we compare it with 
the solution of the system of coupled equations for linear perturba- 
tions <19k numerically integrated with CMBEASY. The comparison 
is shown in Fig. [5] for all the constant coupling models. The accu- 
racy of the linear growth computed from the simulations in fitting 
the theoretical prediction is of the same order for all the values 
of the coupling, and the discrepancy with respect to the numerical 
solution obtained with our modified version of GADGET-2 never 
exceeds a few percent. 



3.3 Our set of N-body simulations 

In our simulations, we are especially interested in the effects that 
the presence of a coupling between DE and CDM induces in the 
properties of collapsed structures, and we would like to understand 
which of these effects are due to linear features of the coupled 
theory, and which due to the modified gravitational interaction in 
the dark sector. This goal turns out to be challenging due to the 
presence of several different sources of changes in the simulation 
outcomes within our set of runs. To summarize this, let us briefly 
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Figure 4. Comparison of the function f(a) with its usual approximation / = f2^ 55 and with the new fit of Eqn. )4U for a ACDM model and for a series of 
coupled DE models. 



discuss in which respect, besides the different gravitational interac- 
tions, the high-resolution simulations listed in Table|2]are different 
from each other: 

• the initial conditions of the simulations are generated using a 
different matter power spectrum for each model, i.e. the influence 
of the coupled DE on the initial power spectrum is taken into ac- 
count and this means that every simulation will have a slightly dif- 
ferent initial power spectrum shape; 

• the amplitude of density fluctuations is normalized at z = 
for all the simulations to as = 0.796, but due to the different shapes 
of the individual power spectra the amplitude of density fluctua- 
tions at the present time will not be the same in all simulations at 
all scales; 

• the initial displacement of particles is computed for each simu- 
lation by scaling down the individual power spectrum amplitudes 
as normalized at z = to the initial redshift of the simulations 
(zi = 60) by using for each simulation the appropriate growth 
function. This results in a lower initial amplitude for more strongly 
coupled models; 



• hydrodynamical forces are acting on baryon particles in all 
the four fully self-consistent simulations (ACDM, RP1, RP2, RP5), 
and therefore differences in the evolution of the dark matter and 
baryon distributions might be due to a superposition of hydrodyna- 
mics and modified gravitational interaction; 

• non-adiabatic processes like e.g. radiative cooling, star forma- 
tion, and feedback are not included in any of the simulations pre- 
sented in this work. 

In order to try to disentangle which of these differences cause 
significant changes in our results, we decided to run three further 
test simulations in which, in turn, some of the new physics has been 
disabled. 

• In the two simulations labelled as "NO-SPH" (ACDM- 
NO-SPH, RP5-NO-SPH), we disabled hydrodynamical SPH 
(Smoothed Particle Hydrodynamics) forces in the code integration. 
We can then compare a ACDM model with a strongly coupled 
model treating both baryons and cold dark matter particles as colli- 
sionless particles. The differences in the dynamics will then be due 
only to the different gravitational interaction implemented in the 
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Figure 5. Evolution of the growth function with redshift for the seven mod- 
els of coupled DE investigated with the low-resolution simulations ran with 
our modified version of GADGET-2. The solid lines are the total growth 
functions as evaluated numerically with CMBEASY, while the triangles are 
the growth function evaluated from the simulations. The relative accuracy 
in reproducing the theoretical prediction is of the order of a few percent 
irrespective of the coupling value j3 c . 



RP5 model. However, the shape and amplitude of the initial power 
spectrum for the two simulations are still different; 

• In the simulation labelled RP5-NO-GF, we ran a RP5 cos- 
mological model using as initial conditions the same file we used 
for the ACDM run. This means that no effect arising in this simu- 
lation compared to ACDM can be due to different initial conditions, 
i.e. due to the differences in the shape and amplitude of the initial 
power spectra that are present in the other simulations. 



4 RESULTS 

We now describe our results for the effect of the coupling between 
DE and CDM on nonlinear structures. As first basic analysis steps 
we apply the Friend s-of-Friends (FoF) and SUBFIND algorithms 
(Sprin gel et alj[200lh to identify groups and gravitationally bound 
subgroups in each of our simulations. Given that the seed used for 
the random realization of the power spectrum in the initial condi- 
tions is the same for all the seven simulations, structures will form 
roughly at the same positions in all simulations, and it is hence pos- 
sible to identify the same objects in all the simulations and to com- 
pare their properties. Nevertheless, due to the different timestep- 
ping induced by the different physics implemented in each run, and 
the slightly different transfer functions, objects in the different sim- 
ulations can be slightly offset from each other. We therefore apply a 
selection criterion and identify objects found in the different simu- 
lations as the same structure only if the most bound particle of each 
of them lies within the virial radius of the corresponding structure 
in the ACDM simulation. If this criterion is not fulfilled for all the 
different simulations we want to compare, we do not consider the 
corresponding halo in any of the comparison analysis described be- 
low. We restrict this matching procedure to the 200 most massive 
halos identified by the FoF algorithm, which have virial masses 
ranging from 4.64 x 1O 12 /i _1 M to 2.83 x 10 14 /i~ 1 A/ Q . Note 
that depending on the specific set of simulations we consider in 



our comparative analysis, this can result in small differences in the 
number of halos included in each of the comparison samples. 



4.1 Halo mass function 

For the four fully self-consistent high-resolution simulations listed 
in Table |2] (ACDM, RP1, RP2, RP5), we have computed the halo 
mass function based on the groups identified with the Friends-of- 
Friends (FoF) algorithm with a linking length of A = 0.2 X d, 
where d is the mean particle spacing. It is important to recall that 
our simulations have the same initial random phases and are nor- 
malized to the same erg today, but the shapes of the input power 
spectra are slightly different for each simulation. In Fig. [6] we plot 
the cumulative mass functions for the four simulations at different 
redshifts. Remarkably, the mass functions of all the cosmological 
models and at all th e different red s hifts c onsidered are well fit by 
the formula given in Ijenkins et al] feOOll) . provided it is evaluated 
with the actual power spectrum and the correct linear growth factor 
for the corresponding cosmological model. The usual mass func- 
tion formalism hence continues to work even for coupled DE cos- 
mologies, a r esult in line with recent findings for early dark energy 
cosmologies dGrossi & Springel 200& lFrancis et al.ll2008dlbh . 

We also plot in Fig. [7] the multiplicity function, defined 
as the derivative of the mass function with respect to the 
mass (M 2 /p ■ dn(< M)/dM), for each simulation at differ- 
ent redshifts. This more sensitive representation of the mass 
fu nction reveals a slightly better agreement with the formula 
by ISheth & Tormenl j 19991) compared with that of Ijenkins et al.l 
which are both overplotted for a direct comparison. 



4.2 Matter power spectrum 

The presence of a long-range fifth-force acting only between CDM 
particles induces a linear bias on all scales between the amplitude 
of density fluctuations in baryons and CDM. These density fluctua- 
tions, in fact, start with the same relative amplitude in the initial 
conditions of all the simulations, and then grow at a different rate 
due to the presence of the extra force. Such a bias is then easily 
distinguishable from the hydrodynamical bias that develops only 
at small scales as structure evolves. This effect is clearly visible 
in our simulations, as it can be seen from the baryon and CDM 
power spectra at z = in the four different cosmologies we ana- 
lyze (Fig. [8}: the density power spectra of baryons and CDM end 
up with a different amplitude on all scales at z = in the coupled 
DE models. CDM always has a larger amplitude, and the difference 
grows with increasing coupling f3 c . 

In order to disentangle this effect from the small scale bias due 
to pressure forces acting on baryons, we can make use of our "NO- 
SPH" simulations to show the result if only the fifth-force effect is 
included. This is shown in the last two panels of Fig. [8] 

In order to make the effect described above even more visible, 
and to better show the difference of the hydrodynamical bias from 
the gravitational one induced by the coupled DE component, we 
also plot in Fig. [9] the ratio of the power A 2 (fc) = P(k)k 3 /2n 
in baryons to that in CDM at different redshifts. In these plots, we 
have corrected all the curves for a spurious effect on small scales 
due to the mass difference between baryon and CDM particles in 
the simulations that induces a small drop in the baryon power. This 
effect is of purely numerical origin and could be easily removed 
by using particles of equal mass for the two cosmological matter 
species in the N-body runs. 
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Figure 6. Cumulative mass functions for the four fully self-consistent high-resolution simulations of the coupled DE models. The four differently-coloured 
solid lines in each figure represent the cumul ative mass function at four different redshifts in each of the investigated models. The dot-dashed lines are 
the corresponding predictions according to the Jenkins et al. 12001) formula, computed for each simulation with the appropriate growth function and power 
spectrum. 



Let us now briefly comment on the plots of Fig. [9] The curves 
represent the bias between the baryon and the CDM density fluctu- 
ation amplitudes as a function of the wave number. A constant bias 
of 1.0 at all redshifts is the expected result for two collisionless 
matter species interacting with the same gravitational strength, and 
this is what we find for our ACDM-NO-SPH simulation (dark blue 
curve). In the ACDM simulation, we notice that this value of 1.0 
is maintained at all redshifts only for large scales, while on smaller 
scales, as structures evolve, the collisional nature of baryon parti- 
cles progressively induces a drop of this ratio (black curve). The 
same behaviour is seen for all the other hydrodynamic simulations 
(RP1, RP2, RP5). However, in the latter cases, the large-scale bias 
is always smaller than 1 .0, already at high redshifts, and it decreases 
with increasing values of the CDM coupling j3 c , as expected. This 
is the gravitational bias that appears also at large scales in Fig. [8] 

Particularly interesting is then again the case of the RP5-NO- 
SPH simulation (dark-green curve), as compared to the ACDM- 
NO-SPH one, because it allows us to disentangle the hydrodynamic 
effects from the effects due to the coupled DE extra physics. In this 
case we find that the bias of RP5-NO-SPH perfectly overlaps with 



the one of the other two RP5 simulations at high redshifts, while 
at lower and lower redshifts, the small scale behaviour is progres- 
sively more and more different: as expected the absence of hydro- 
dynamic forces acting on baryon particles induces a larger value 
of the bias, which is now solely due to the different gravitational 
strength felt by the two particle species. It is however very inter- 
esting to notice that the bias does not keep the large-scale linear 
value at all scales, as it is the case for the ACDM-NO-SPH run, 
but evolves towards lower and lower values for smaller and smaller 
scales. This clearly shows that nonlinearities enhance the effect of 
the coupling on the growth of overdensities in the two differently 
interacting matter species. 



4.3 Halo density profiles 

Applying the selection criterion described above to our four self- 
consistent simulations (ACDM, RP1, RP2, RP5) we select among 
the 200 most massive groups identified at z — for each run 74 
objects that can be considered with certainty to be the same struc- 
ture in the different simulations. For these 74 halos we compute 
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Figure 7. Multiplicity functions for the four high-resolution simulations of the interacting DE models. The four differently-coloured sets of data points are the 
multiplicity functions evaluated i n equally spaced loga rith mic mass bins at four diff erent redshifts. The dot-dashed and dashed lines represent the predictions 
for the multiplicity function from Jenkins et al. (2001) and Sheth & Tormen 1 1999), respectively, computed for eac h simulation with the ap propriate growth 
function and power spectrum. The comparison clearly shows a slightly better agreement with the fitting function bv lSheth & Tormenl fl999T) . in particular at 
high redshift. 



the spherically averaged density profiles of CDM and baryons as a 
function of radius around the position of the particle with the min- 
imum gravitational potential. 

Interestingly, the halos formed in the coupled DE cosmolo- 
gies show systematically a lower inner overdensity with respect 
to ACDM, and this effect grows with increasing coupling. This is 
clearly visible in Fig. [10] where we show the density profiles of 
CDM and baryons in the four different cosmologies for a few se- 
lected halos of different virial mass in our sample. We remark that 
this result is clearly incompatibl e with the es s ential ly opposite be- 
haviour previously reported by Macc ubet al.l ^20041) . and deserves 
a more detailed discu ssion. 

Unlike found in iMaccio et all {2004), all the 74 halos in our 
comparison sample ha ve density profiles t hat are well fitted by the 
NFW fitting function dNavarro et al]|l997h 

8^1 = *! , (42) 

Pcrit (r/r s )(l + r/r s ) 2 ' 

independent of the value of the coupling. Here S* is a parame- 
ter that sets the characteristic halo density contrast relative to the 



critical density p cr it- The scale radius r s increases for each halo 
with increasing coupling /3 C , and becomes larger than that found in 
ACDM. In other words, the halos become less concentrated with 
increasing coupling. 

For example, for the four halos shown in Fig.[l0] the scale ra- 
dius grows with increasing coupling by roughly 10% to 35% going 
from ACDM to RP5, as listed in Table g] and shown in Fig. [TT] 
The effect is not strong enough to fully address the "cusp-core" 
problem, but clearly shows how the interaction between DE and 
CDM can produce shallower halo density profiles with respect to a 
ACDM cosmology with the same cosmological parameters. This 
result goes then in the direction of alleviating the problem, and 
therefore opens up new room for the study of dark interactions 
as opposite to previous claims. In fact, although the models pre- 
sented here span over the full obs ervationally allow ed parameter 
spa ce for constant-coup ling models JBean et al. ( 2008), but see also 
e.g. lLa Vacca et al.l |2009)), it is well conceivable that more realis- 
tic scenarios with a variable coupling strength could produce larger 
effects in the nonlinear regime without running into conflict with 
present observational bounds from linear probes. We defer to fu- 
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Figure 8. Power spectra of CDM (black line) and baryons (blue line) at z = for the set of coupled DE models under investigation. The appearance of a 
bias between the two distributions, which grows with increasing coupling /3 C , is clearly visible at the large scale end of the plots. The last two panels show the 
comparison of a ACDM and a coupled DE cosmology with /3 C = 0.2 in absence of hydrodynamic forces acting on baryons. In these two panels, the bias on 
all scales is purely due to the interaction of CDM with the DE scalar field (f>- 
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Figure 9. Ratio of the power spectra of baryons and CDM as a function of wavenumber for the set of high-resolution simulations ran with our modified version 
of GADGET-2, for five different redshifts. The linear large-scale bias appears already at high redshifts, while at lower redshifts the hydrodynamic forces start 
to suppress power in the baryon component at small scales. In absence of such hydrodynamic forces the progressive enhancement of the large scale bias at 
small scales for the RP5-NO-SPH run (light green curve) as compared to the completely flat behaviour of the ACDM-NO-SPH simulation (blue curve) - 
where no bias is expected - shows clearly that nonlinearities must increase the effect of the coupling on the different clustering rates of the two species. All 
the curves have been corrected for a spurious numerical drop of the baryonic power at small scales as described in the text. 
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Figure 10. Density profiles of CDM (solid lines) and baryons (dot-dashed lines) for four halos of different mass in the simulation box at z = 0. The vertical 
dot-dashed line indicates the location of the virial radius for the ACDM halo. The decrease of the inner overdensity of the profiles with increasing coupling is 
clearly visible in all the four plots. 



ture work the investigation of such models, but we stress here that 
our present results constitute the first evidence that the nonlinear 
dynamics of generalized coupled cosmologies might provide a so- 
lution to the "cusp-core" problem. 



4.4 Halo concentrations 

For all the 200 most massive halos found in each of our four fully 
self consistent simulations we compute halo concentrations as 



rgoo 
r 3 



(43) 



based on our NFW fits to the halo density profiles. Here r2oo is the 
radius enclosing a mean overdensity 200 times the critical density. 
Note that here no further selection criterion is applied, and the con- 
centration is computed for all the 200 most massive halos in each 
simulation. 

Consistently with the trend found for the inner overdensity in 
the halo density profiles and for the evolution of the scale radius 
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Figure 11. Relative evolution with respect to ACDM of the scale radius r s 
for the four halos plotted in Fig.[l0]as a function of coupling f) c . 
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Table 4. Evolution of the scale radius r s for the four halos shown in Fig.[lO]with respect to the corresponding ACDM value. The trend is towards larger values 
of r B with increasing coupling /3 C , with a relative growth of up to 36% for the largest coupling value /3 C = 0.2. 



with coupling, we find that halo concentrations are on average sig- 
nificantly lower for coupled DE models with respect to ACDM, 
and the effect again increases with increasing coupling f3 c . This 
behaviour is shown explicitly in the left panel of Fig. [T2j where 
we plot halo concentrations as a function of the halo virial mass 
M200 for a series of our high-resolution simulations. In the stan- 
dard interpretation, the halo concentrations are thought to reflect 
the cosmic matter density at the time of formation of the halo, lead- 
ing to the association of a larger value of the concentration with 
an earlier formation epoch, and vice versa. In the context of this 
standard picture, the effect we found for the concentrations could 
be interpreted as a sign of a later formation time of massive ha- 
los in the coupled DE models as compared to the ACDM model. 
Such a later formation time could be possibly due to the fact that 
matter density fluctuations start with a lower amplitude in the ini- 
tial conditions of the coupled cosmologies with respect to ACDM, 
and this would make them forming massive structures later, despite 
their faster linear growth (as shown in Fig. [5}- 

However, we can demonstrate that this is not the case, just 
making use of our RP5-NO-GF simulation, in which the Universe 
evolves according to the same physics as RP5, but starting with the 
identical initial conditions as used for the ACDM run. Therefore, 
any difference between these two simulations can not be due to the 
initial amplitude of fluctuations. The evolution of halo concentra- 
tions with mass for this run is also plotted in Fig. [T2] (dark blue 
curve), and shows a very similar behaviour to the RP5 curve. 

As a cross check of this result, we have repeated the same anal- 
ysis by computing halo concentrations with an independent method 
that circumvents the profile fitting. The concentration can be re- 
lated to two other basic structural properties of a halo, namely its 
maximum rotational velocity Knax, and th e radius at wh i ch thi s 
velocity peak is located, r max - According to lSpringel et alj |2008), 
the concentration can then be related to these two quantities by the 
relation: 



Mean formation redshift for 
different coupled dark energy models 



200 



= 7.213 <5 V 



3 ln(l + c) - c/(l + c) 
where 5v is a simple function of V x 

Vt 



and r n 



8 V = 2 



/ Vmax V 
I /Jo r max/ 



(44) 



(45) 



v Hor max , 

We denote the concentrations evaluated in this way as c*, and in- 
clude our results as a function of halo mass in the right panel of 
Fig. [T2] Although not identical in detail, as expected due to the 
different methods used to measure concentrations, the two plots 
of Fig. [12] show the same trend for the evolution of halo concen- 
trations with coupling, and the same independence of this effect 
from the initial conditions of the simulations. The slight down-turn 
of the concentrations curve at the low-mass end that is visible in 
both the two panels of Fig. [12] is an expected effect due to nu- 
merical resolution, as was discussed in detail in iNeto et alj 120071) 
and in ISpringel et all {2008). Although the concentration decrease 
is found to have the same trend and amplitude for the poorly re- 
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Figure 13. Evolution of the mean halo formation redshift 2 * as a function 
of halo mass for the 200 most massive halos in our simulations and for the 
different cosmological models under investigation. The formation redshift 
Zf is defined as the redshift at which the main progenitor of the halo has a 
virial mass equal to half the virial mass of the halo at z = 0. The halos have 
been binned by mass, and the mean formation redshift in each bin is plotted 
as a filled circle. The coloured dashed lines indicate for each simulation the 
spread of 68% of the halos in each mass bin. The highest mass bin is not 
plotted because of its too low number of halos. 



solved low-mass halos as for the larger structures for which no res- 
olution problem is expected, one should keep in mind as an element 
of caution that higher resolution simulations would be required for 
a definitive confirmation of this effect in the low-mass range of our 
comparison sample, i.e. for masses below 10 13 Mq/Ii. 

In order to directly verify that the lower concentrations cannot 
be a consequence of a later formation time we have also computed 
the average formation redshift of the halos in our sample for all the 
four self-consistent simulations, by building merger trees for all the 
halos in our sample, and by following backwards in time the main 
progenitor of each halo until the redshift at which its virial mass is 
only half of the final virial mass of the halo at z = 0. We define the 
corresponding time as the formation redshift Zf of the halo. 

In Figure [T3] we show the evolution of zj as a function of 
halo mass for all our cosmological models. It is evident that mas- 
sive halos in the different cosmologies form approximately at the 
same time, with a slightly earlier formation for the RP5 cosmol- 
ogy, which one might have expected to translate into slightly larger 
values of the concentrations. Therefore we conclude that the un- 
ambiguous trend of lower halo concentrations for larger coupling 
values must be a peculiar feature that arises from the extra physics 
that characterizes the coupled DE cosmologies. A more detailed in- 
vestigation of how this peculiar behaviour arises is hence required 
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Figure 12. Evolution of the mean halo concentration as a function of mass for the 200 most massive halos in our simulations and for the different cosmological 
models under investigation. The conc entrations have been c omputed by directly fitting the halo density profile of each halo with an NFW model (left panel) 
or by using the method introduced by Springel et al. (2008) and described in Eqs. j44!45> (right panel). The halos have been binned by mass, and the mean 
concentration in each bin is plotted as a filled circle. The coloured dashed lines indicate for each simulation the spread of 68% of the halos in each mass bin. 
The highest mass bin is not plotted because of its very low number of halos. The decrease of the mean concentration with increasing coupling appears in the 
same way in both plots. 



in order to understand this phenomenology of the dynamics in cou- 
pled DE cosmologies. 

We perform such an investigation by switching off individu- 
ally the two main effects which could be responsible of the con- 
centrations drop, which are the variation of particle mass and the 
extra velocity-dependent term. To save computational time, we do 
this only for the most strongly coupled model RP5 and only for 
the late stages of cosmic evolution. More specifically, we take one 
of our RP5 simulation snapshots at a given redshift z* , and use it 
as initial conditions file for a new run starting at z = z* down to 
z = in which one of these two effects is switched off. We la- 
bel these simulations as "RP5-NO-MASS" and "RP5-NO-FRIC" 
for the cases where the mass decrease or the velocity-dependent 
term are dropped, respectively. We set z* = 1.5 as a conservative 
choice based on the consideration that, according to our definition 
of formation redshift of a halo, all the halos in our sample have a 
formation redshift z < z* , as shown in Fig. 1 131 

By switching off the mass variation for z < 2*, we find that 
the halo concentrations at 2 = show a slight increase over the 
whole mass range of the sample with respect to the fully self- 
consistent RP5 simulation. This effect is shown in the left panel of 
Fig-UU We interpret this as a sign of the fact that the mass decrease 
reduces the total gravitational potential energy of halos, resulting in 
a modification of their virial equilibrium configuration. 

In fact, if the potential well of a halo gets shallower as time 
goes by as a consequence of the decrease of the mass of its CDM 
content, the system will find itself with an excess of kinetic en- 
ergy, and will therefore expand in order to restore virial equilib- 
rium. Such an expansion is expected to cause a drop of the halo 
concentrations, which we confirm here because switching off this 
mechanisms yields consistently higher concentrations at z = 0. 
However, it is clear from Fig. [14] that this mechanism cannot ac- 
count for the total effect of concentration decrease, but only for a 
small fraction of it. 

We therefore now investigate the other possible origin of this 
effect, i.e. the impact of the velocity-dependent term l |29t on the dy- 



namics of CDM particles. To this end we switch off for 2 < z* the 
additional acceleration arising from the velocity-dependent term 
for coupled particles described by Eqn. J3. l-3t . The outcome of this 
test simulation is shown in the right panel of Fig. [14] the increase 
of the concentrations with respect to the fully self-consistent RP5 
simulation is now much more substantial than in the case of RP5- 
NO-MASS, and shows that the velocity-dependent term is actually 
the dominant mechanism in determining the decrease of halo con- 
centrations and the decrease of the inner overdensity of CDM halos 
discussed above. The interpretation of this effect seems quite un- 
ambiguous: the velocity-dependent term induces an extra accelera- 
tion on coupled particles in the direction of their velocity, and this 
produces an increase of the kinetic energy of the particles, moving 
the system out of its virial equilibrium configuration. The system 
responds by a small expansion and a lowering of the concentration. 

As a further check of this interpretation of our results, we also 
test directly the dynamic evolution of halos to check whether they 
really slightly expand in the presence of DE-CDM coupling. To this 
end, we compute for all the halos in our sample the time evolution 
of the mass and the number of particles contained in a sphere of 
physical radius r — 20 /i _1 kpc centred on the potential minimum 
of each halo. This sphere represents the very innermost part of all 
the halos in our sample at any redshift between 2* and 0, and we 
refer to it as the halo "core"; its mass content is expected to be 
roughly constant for ACDM cosmologies at low redshifts. Indeed, 
we can recover this behaviour for our ACDM simulation by aver- 
aging the evolution of core masses and particle numbers over the 
whole halo sample. On the other hand, for RP5, as expected ac- 
cording to our interpretation, both the mass and the number of par- 
ticles in the halo cores strongly decrease with time. This is shown 
in Fig.Q3] where the solid lines represent the average evolution of 
mass, and the dashed lines represent the average evolution of the 
particle number. Evidently, for ACDM the two curves coincide be- 
cause the mass of the particles is constant and any change of the 
enclosed mass in the core must be due to a change of the number 
of enclosed particles. On the other hand, for RP5, the mass and 
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Figure 14. Evolution of halo concentrations for the same models and the same halo sample as in Fig. 1121 and for an additional test simulation in each of the 
two panels. In the left panel, the simulation RP5-NO-MASS shows the effect of switching off the mass correction for z < z* ~ 1.5: there is a small but 
systematic increase of average halo concentrations over the whole mass range. In the right panel, the simulation RP5-NO-FRIC shows the effect of switching 
off in the same redshift interval the velocity-dependent term. The increase of concentrations is in this case much more consistent and accounts for a large 
fraction of the total concentration reduction of RP5. 



the particle number behave differently due to the mass variation of 
CDM particles. The decrease of the number of particles contained 
in the core can be interpreted as a manifestation of an expansion 
of the halos. Moreover, if we compute the same evolution for our 
RP5-NO-FRIC simulation, we find an almost constant evolution 
of the core particle number and a very weak decrease of the core 
mass due to the variation of CDM particle mass. This result also 
confirms our interpretation concerning the origin of the decrease 
of concentrations: the velocity-dependent term is the most relevant 
mechanism for inducing halos expansion at low redshifts, and as a 
consequence the decrease of the inner overdensity of CDM halos 
and of their concentration. 

While further investigations of these effects are certainly re- 
quired in order to understand all the potential phenomenological 
features of interacting DE models, our conclusion that a coupling 
between DE and CDM produces less peaked halo density profiles 
and lower halo concentrations seems to be quite robust based on 
the analysis of our simulations that we have discussed here. We 
note again that ou r findings are i n star k contrast with the results of 
previous work bv lMaccio et al. (2004) who found for coupled DE 
models a strong increase in concentration and density profiles in the 
centre that more steeply rise than ACDM. The effects we find go in 
the direction of less "cuspyness" of halo density profiles, which is 
preferred by observations and thus in fact opens up new room for 
the phenomenology of interacting DE models. 



4.5 Integrated bias and halo baryon fraction 

The extra force felt by CDM particles induces, as we have already 
seen for the evolution of the matter power spectrum, a bi as in the 
evolution of density fluctuations of baryons and CDM |Mainini 
120051: iMainini & Bonome"tt3l200d : lManera & Motall2006l) . We can 
then use our selected halo sample to test the evolution of this bias 
from the linear regime already probed by the power spectrum on 
large scales to the highly nonlinear regime in the centre of massive 
collapsed structures. We then test the evolution of the integrated 
bias 
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Figure 15. Evolution with respect to the scale factor a of the average mass 
(solid) and of the average number of particles (dashed) enclosed in a sphere 
of physical radius r = 20h~ 1 kpc centred on the potential minimum of 
each halo in our sample. The curves are normalized at a = 0.48 (z ~ 1) 
and show the expected flat behaviour for the ACDM case (black line) for 
which the solid and the dashed curves coincide due to the constancy of the 
mass of particles. For the RP5 case (light blue curves), there is a strong 
decrease in time of both mass and particle number, which clearly illustrates 
the expansion of RP5 halos with respect to the ACDM case. By switching 
off the extra velocity-dependent force acting on CDM particles (RP5-NO- 
FRIC, light green curves), an almost flat behaviour is recovered again for 
the particle number, while the decrease of mass is now due to the particle 
mass variation - which is still in place for this simulation - on top of the 
particle number evolution. This plot therefore clearly shows that the extra 
physics of coupled DE cosmologies induces an overall expansion of CDM 
halos at low redshifts, and clearly identifies in the velocity-dependent 

term the leading mechanism that produces this expansion. 
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B(< 



Pb(< r) - p b 



Pc ' 



(46) 



pb Pc(< r) 

as defined in IMaccio et al] d2004l) . where pb(< r) and p c (< r) 
are the densities within a sphere of radius r around the potential 
minimum of a halo, for baryons and CDM, respectively. Following 
IMaccio et alj d2004l) . we have not used the innermost part of the ha- 
los (r < 10/i _1 kpc ~ 3 x e s ) in order to avoid potential resolution 
problems. 

In Figure [16] we show the evolution of the bias for four se- 
lected halos of our sample with similar masses as the ones shown 
in Fig. [10] It clearly appears that the bias is considerably enhanced 
in the nonlinear regime, while at large scales it converges to the 
linear value evaluated from the power spectrum amplitude on large 
scales, represented in Fig.QJojby the horizontal dashed lines. 

However, also in this case this effect could be due only to the 
presence of hydrodynamical forces acting on the baryons, and may 
not really be caused by the fifth-forces from the coupled DE scalar 
field, as we can infer from the fact that also the ACDM curve, 
where no coupled DE is present, shows a departure from the large 
scale value of 1.0 when approaching the centre of the halos. Once 
again we make use of our additional test simulations ACDM-NO- 
SPH and RP5-NO-SPH in order to disentangle the two effects. In 
Fig- El] we show the same four plots as in Fig.[T6]for the two simu- 
lations without hydrodynamic forces, and the appearance of a non- 
linear bias imprinted only by the coupled DE scalar field acting on 
CDM particles is then absolutely evident. On the other hand, the ab- 
sence of any bias, as expected, in the ACDM-NO-SPH run shows 
clearly that no major numerical problems can be responsible for the 
effect in the RP5-NO-SPH simulations. This evolution of the bias 
B(< r) with radius for massive halos is in good agreement with 
that found by |Macci6et~al]d2004l) despite the starkly different be- 
havior of the individual overdensities in baryons and CDM found 
in the two works. 

It is interesting that the above effect produces a baryon deficit 
in virialized halos, i.e. they contain fewer baryons than expected 
based on their mass and the universal cosmological baryon fraction. 
In particular, this means that one can not expect that baryon frac- 
tions determined through X-ray measurements in clusters would 
yield the cosmological value. It is important to stress again here 
that our approach does not include non-adiabatic processes (like 
e.g. star formation, cooling, and possible feedback mechanisms) 
to the physics of the baryonic component. Therefore our outcomes 
cannot be expected to provide a prediction of baryon fractions to be 
directly confronted with observations. However, this approach has 
the advantage to clearly show in which direction the DE coupling 
affects the baryon budget of collapsed structures which would then 
be available for such radiative processes. In order to give a rough 
estimate of the magnitude of this effect we compute the baryon 
fraction within the virial radius r2oo of all the halos in our sample 
defined as 



h 



M b (< rgoo) 
A/tot (< r 20 o) 



(47) 



for our four fully self-consistent simulations. We plot in Fig.l 1 8las 
a function of halo virial mass the relative baryon fraction defined 



fb 



(48) 



For the ACDM case, our results for the evolution of Yb are consis- 
tent with the value of Yj, ~ 0.92 found by t he Santa Barbara Clus- 
ter Comparison Project dFrenk et al.|[l999l) , and with the more re- 
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Figure 18. Evolution with virial mass M200 of the relative baryon fraction 
Yf, within the virial radius r200 of all the halos in our sample. The coloured 
diamonds represent the relative baryon fraction of each single halo, while 
the filled circles and the coloured curves show the behaviour of the mean 
relative baryon fraction in each mass bin for the four fully self-consistent 
high-resolution simulations. A decrease of Yj, with increasing coupling is 
clearly visible both in the distribution of the individual halos and in the 
averaged curves. 



cent results of lEttori et ail d2006h and lGottloeber & Yepesl d2007l) . 
while for the coupled models the relative baryon fraction shows a 
progressive decrease with increasing coupling, down to a value of 
Y b ~ 0.86 - 0.87 for the RP5 case. 

It is also important to notice that this effect is always to- 
wards lower baryon fractions in clusters with respect to the cos- 
mological value. This could in fact alleviate tensions between the 
high baryon abundance estimated from CMB observations, and 
the somewhat lower value s inferred from detailed X-ray obser- 
vations of galaxy clusters dVikhlinin et al J 2006; Mc Carthy et al.l 
120071; iLaRoaue et ai]2006l : lAfshordi et alj2007l) . 



5 CONCLUSIONS 

We have investigated coupled DE cosmologies, both with respect 
to their expected background and linear perturbation evolution, as 
well as in their predictions for the nonlinear regime of structure 
formation. To do so we have developed and tested in this work a 
modified version of the cosmological N-body code GADGET-2 sui- 
table for evolving these kinds of cosmological models. The nume- 
rical implementation we have developed is in fact quite general and 
not restricted to the simple specific models of coupled quintessence 
that we have investigated in this paper. Instead, it should be well 
suited for a much wider range of DE models. We also note that 
the ability to selectively enable or disable each of the modifications 
discussed above, makes the code suitable for cosmological models 
that are unrelated to the coupled DE scenario but require similar 
new degrees of freedom that our implementation allows. These are: 

(i) the expansion history of the Universe can be specified ac- 
cording to any desired evolution of the Hubble rate as a function of 
the scale factor a; 

(ii) a global variation in time of the gravitational constant and/or 
a variation in time of the gravitational strength for each individual 
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Figure 16. Evolution of the integrated bias B(< r) for the four fully self-consistent high-resolution simulations and for four selected halos of different mass 
in our sample. The horizontal dashed lines indicate the value of the large scale linear bias as evaluated from the power spectrum amplitudes of baryons and 
CDM. The vertical black dot-dashed line shows the position of the virial radius for the ACDM halo in the sample. The drop of the value of B(< r) in the 
innermost regions of the halos is evident but in these runs is given by a superposition of effects due to hydrodynamical forces and to the modified gravitational 
interaction. On large scales, the bias tends to converge to the linear value, as expected. 



matter species. This includes the possibility to have a long range 
repulsive interaction between different particle species; 

(iii) variation in time of the particle mass of each individual mat- 
ter species; 

(iv) extra velocity-dependent terms in the equation of motion for 
each individual matter species. 

With this implementation we have investigated the effects of 
coupled DE models with a constant coupling /3 C to the CDM fluid 
on structure formation. We have shown that the halo mass func- 
tion is modified in coupled DE models, bu t can still be well fi tted 
at different redshifts by the Jenkins et al. ( Jenki ns et alj|200lh fit - 
ting formula, or by the Sheth & Tormen fsheth & Tormen Jl999h 
formula, which yields a moderately better agreement, especially at 
2 > 0. 

We have confirmed the analytic prediction that density fluc- 
tuations in baryons and CDM will develop a bias on all scales due 
to the presence of a fifth-force acting only between CDM particles. 
We have also shown that in addition to this the bias is enhanced 
when moving from the linear regime of very large scales to smaller 
and progressively more non-linear scales. 

We have investigated the evolution of the bias between 
baryons and CDM overdensities down to the very nonlinear regime 
found in the inner part of collapsed objects, in the same fashion 



as described in lMaccio et alj ^20041) . We found here similar results 
with this previous work, namely an enhancement of the bias in the 
nonlinear region within and around massive halos. We also recover 
from this analysis the large scale value of the linear bias computed 
from the power spectrum when integrating the bias function up to 
very large radii from the centre of CDM halos. The enhancement of 
the bias in highly nonlinear structures has an impact on the deter- 
mination of the baryon fraction from cluster measurements, and we 
have computed for all our halos, without including non-adiabatic 
processes, the evolution of this fraction within the virial radius r2oo 
with coupling, finding that the baryon fraction is reduced with in- 
creasing coupling by up to ~ 8 — 10% with respect to ACDM for 
the largest coupling value. 

We have also investigated the effect of the coupling on the halo 
density profiles. We find that they are remarkably well fit over the 
resolved range by the NFW formula for any value of the coupling. 
There is a clear trend of a decrease of the inner halo overdensity 
with respect to ACDM with increasing coupling (or, equivalently, 
an increase of the scale radius r s for increasing coupling). This 
result confli cts with previous cl aims for the same class of coupled 
DE models jMaccio et al]2004h . 

Using a number of special test simulations, we have identified 
the origin of this effect of reduced halo concentrations for increas- 
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Figure 17. Evolution of the integrated bias B(< r) for the two high-resolution simulations without hydrodynamical forces on baryon particles for the same 
four halos shown in Fig. 1 161 The enhancement of the bias due to the extra scalar force in the core of highly nonlinear structures appears here clearly. 



ing coupling. It actually arises from a combination of two peculiar 
features that the coupling introduces in the Newtonian limit of grav- 
itational dynamics. The first of these is the decrease of CDM parti- 
cle mass with time, which causes the total potential energy of a halo 
to decrease, and hence effectively moves the system to a configu- 
ration where an excess of kinetic energy is present relative to virial 
equilibrium. The second one is the additional velocity-dependent 
term, which directly raises the total kinetic energy density of halos 
by accelerating their particles in the direction of their peculiar ve- 
locity. Both of these effects cause a halo to slightly expand in order 
to restore virial equilibrium, and this reduces the halo concentra- 
tion. 



In conclusion, we have developed a general numerical im- 
plementation of coupled dark energy models in the GADGET-2 
code. We have then performed the first fully self-consistent high- 
resolution hydrodynamic N-body simulations of interacting dark 
energy models with constant coupling, and carried out a basic 
analysis of the non-linear structures that formed. Interestingly, we 
found that a larger coupling leads to a lower average halo concen- 
tration. Furthermore, both the baryon fraction in massive halos and 
the inner overdensity of CDM halos decrease with increasing cou- 
pling. These effects alleviate the present tensions between observa- 
tions and the ACDM model on small scales, implying that the cou- 
pled DE models are viable alternatives to the cosmological constant 
included in standard ACDM. 
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